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What are the relevant timescales of neural encoding in the brain? This question is 
commonly investigated with respect to well-defined stimuli or actions. However, neurons 
often encode multiple signals, including hidden or internal, which are not experimentally 
controlled, and thus excluded from such analysis. Here we consider all rate modulations 
as the signal, and define the rate-modulations signal-to-noise ratio (RM-SA/ff) as the 
ratio between the variance of the rate and the variance of the neuronal noise. As the 
bin-width increases, RM-SA/f? increases while the update rate decreases. This tradeoff 
is captured by the ratio of RM-SA/ff to bin-width, and its variations with the bin-width 
reveal the timescales of neural activity. Theoretical analysis and simulations elucidate 
how the interactions between the recovery properties of the unit and the spectral 
content of the encoded signals shape this ratio and determine the timescales of neural 
coding. The resulting signal-independent timescale analysis (SITA) is applied to investigate 
timescales of neural activity recorded from the motor cortex of monkeys during: (i) 
reaching experiments with Brain-Machine Interface (BMI), and (ii) locomotion experiments 
at different speeds. Interestingly, the timescales during BMI experiments did not change 
significantly with the control mode or training. During locomotion, the analysis identified 
units whose timescale varied consistently with the experimentally controlled speed of 
walking, though the specific timescale reflected also the recovery properties of the unit. 
Thus, the proposed method, SITA, characterizes the timescales of neural encoding and 
how they are affected by the motor task, while accounting for all rate modulations. 

Keywords: timescales, rate modulation, doubly stochastic Poisson processes, rate coding, brain-machine interface, 
signal to noise ratio 



INTRODUCTION 

The timescales of neural encoding is important for understand- 
ing neural processing and for successfully interpreting recorded 
neural activity (Shadlen and Newsome, 1994; de Ruyter van 
Steveninck et al, 1997; Borst and Theunissen, 1999; Jacobs et al, 
2009). In the context of rate-coding, we focus on characteriz- 
ing those timescales, independent of which signals modulate the 
rate. Those may include: (i) encoded signals under experimen- 
tal control, (ii) encoded task-relevant signals that are not under 
experimental control and thus may vary from trial to trial, and 
(iii) task-irrelevant signals. The last two groups are referred to 
as "hidden" signals, and we are especially interested in the sec- 
ond group of hidden task-relevant signals. Those may include, 
for example, the estimated state or estimation error during reach- 
ing movements (Desmurget and Grafton, 2000; Wolpert and 
Ghahramani, 2000; Krigolson and Holroyd, 2006; Shadmehr and 
Krakauer, 2008). 

The potential contribution of hidden signals to trial-to-trial 
variability in neuronal responses has been noted in the con- 
text of different tasks, including skilled reaching movements 
(Churchland et al., 2006a,b; Mandelblat-Cerf et al, 2009) motor 
adaptation (Mandelblat-Cerf et al., 2009) and decision making 



(Churchland et al, 2011). During skilled reaching movements, 
trial-to-trial variations in preparatory neural activity were pre- 
dictive of variations in reach, and were interpreted to reflect 
variations in movement planning (Churchland et al, 2006a). 
During novel visuo-motor tasks, trial-to-trial variability in spik- 
ing activity increased in early adaptation stages before return- 
ing to initial levels at the end of learning, and was interpreted 
to reflect enhanced exploration or adjustments of the internal 
models (Mandelblat-Cerf et al., 2009). During decision mak- 
ing, trial-to-trial variations in the underlying rate were shown 
to increase until the decision was made, in agreement with the 
hypothesis that they encode trial-specific evidence accumulation 
(Churchland etal, 2011). 

Ideally we should quantif)? and analyze the contributions of the 
first two groups of signals, i.e., only those that are task-relevant, 
whether under experimental control or hidden. However, this 
is hampered by the inability to control or measure hidden 
task-relevant signals. Current methods for analyzing neuronal 
timescales, which have focused mainly on sensory neurons, are 
based on averaging the response over repeated trials with the 
same stimuli (Borst and Theunissen, 1999; Warzech and Egelhaaf, 
1999; Butts et al., 2007), or computing the mutual information 
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between known stimuli and their reconstruction (Rieke et al., 
1997; Dimitrov and Miller, 2000). In either case, the analysis 
accounts only for signals that are under experimental control, i.e., 
the first group of signals. 

Neurons in the motor system are known to encode a large 
number of signals (Georgopoulos, 2000; lohnson et al, 2001; 
Scott, 2003), and are also expected to encode hidden signals, such 
as the estimated state and errors (Desmurget and Grafton, 2000; 
Wolpert and Ghahramani, 2000; Krigolson and Holroyd, 2006; 
Shadmehr and Krakauer, 2008). Timescale analysis that is based 
on synchronized averaging would miss the contribution of those 
signals that are not under experimental control, including hidden 
task-relevant signals. This is illustrated in Figure lA (detailed in 
Section Rate-modulations Signal-to-noise ratio of DSPP), where 
the rate is assumed to encode both a movement related sig- 
nal and a hidden signal: synchronized averaging captures the 
rate-encoded movement but not the rate-encoded hidden sig- 
nal. Instead, we focus on analyzing the timescales associated with 
the underlying rate modulations, independent of which signals 
modulate the rate. The proposed method considers asynchronous 
sequences of random reaching movements (Figure IB, upper 
panel) so the resulting rate (Figure IB, lower panel) is a stationary 
process whose variance (dashed black) captures the variance of 
both of the encoded signals (dashed black lines. Figure lA upper 
panel and Figure IB lower panel, respectively). Thus, in contrast 
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FIGURE 1 I Illustration of the averaging approach (A) vs. the 
proposed approach (B) when the rate encodes both the movement 
and a hidden signal. Averaging methods are based on repeating tine 
same movement (A, upper panel) and synchronizing the resulting neura 
activity and hence the rate (A, lower panel). Synchronized average (solid 
black) is considered the signal, while the variance of the hidden signal 
(dashed black) is considered noise. The proposed method considers 
asynchronous sequences of random reaching movements (B, upper 



with standard averaging methods, the proposed method captures 
also the contribution of hidden signals. Admittedly, the hidden 
signals may include not only task- relevant (group 2) but also task- 
irrelevant (group 3) hidden signals. Hence, the method is best 
used in comparing the timescales across different task conditions. 
Here we demonstrate the power of this method by compar- 
ing the timescales in neural activity during different phases of 
experiments with brain-machine interfaces (BMIs) and during 
locomotion at different speeds. 

Timescales are assessed under the sole assumption that spike 
trains are realizations of dead-time modified doubly stochastic 
Poisson processes (DSPP) (Snyder, 1975; Cox and Isham, 1980; 
Johnson, 1996; Gabbiani and Koch, 1998). DSPPs are the simplest 
point processes that can describe rate modulations by stochas- 
tic signals. This model is pertinent for describing spike trains 
recorded during reaching and walking when the rate might be 
modulated by a number of biologically relevant stochastic signals. 
These signals may include not only measurable variables, such 
as the velocity of movement, but also internal or hidden signals 
(e.g., internal state estimations; Wolpert and Ghahramani, 2000; 
Shadmehr and Krakauer, 2008), which cannot be measured or 
controlled directly (Zacksenhouse et al., 2007). The effect of dead- 
time is explicitly investigated to account for absolute refractory 
period. The effects of other deviations from the DSPP assumption 
are demonstrated via simulations. 



B Asynchronized random reach 
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panel) so the resulting rate (B, lower panel) is a stationary process 
whose variance (dashed black) captures the variance of both of the 
encoded signals Idashed black lines, (A) upper panel and (B) lower 
panel, respectively]. Each movement has a minimum jerk profile over 
500 ms. Random reaching movements (B) are 2 s long sequences of 
reaching movements between random targets. Colored traces: 
representative samples. Black solid and dashed lines: ensemble mean 
and variance computed across an ensemble of 500 samples. 
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Within the DSPP model, the modulated rate is considered 
as the signal conveyed by the unit. Only rate-conditioned spike 
count variations are considered noise (Churchland et al., 2011). 
Thus, we define the rate-modulations signal-to-noise ratio (RM- 
SNR) as the ratio between the variance of the rate and the 
variance of the Poisson noise (Section Rate-modulations Signal- 
to-noise ratio of DSPP). Theoretical analysis, presented in Section 
Bandwidth Effect and Appendix A, indicates that at short bin- 
widths RM-SM^ should increase linearly with the bin-width, 
while at long bin-widths it should saturate. The analysis is 
extended in Section Dead-time Effect and Appendix B to include 
the effect of refractory period (dead-time modified DSPP). 
Section Simulations describes the simulations that are used to 
demonstrate the proposed methods. Section Timescales Analysis 
proposes to capture the trade-off between RM-SNR and update- 
rate by the ratio between RM-SNR and bin-width. The variations 
of this ratio with the bin-width are related to the interaction 
between the refractory properties of the unit and the spectral 
properties of the encoded signals, and the bin-widths at which 
this ratio peaks are related to the timescales of neural encoding. 
The effects of deviations from the DSPP assumption are evaluated 
by analyzing simulated realizations of doubly stochastic Gamma 
processes (DSGPs). 

After demonstrating the analysis on simulated realizations of 
both DSPPs and doubly stochastic Gamma processes (DSGPs), 
it is applied to investigate the timescales of the neural activity 
recorded from cortical motor units during two experiments with 
Monkeys (described briefly in Section Experimental Methods). 
Sections SNR During BMI Experiments and Timescale During 
BMI Experiments present the analysis of spike-trains recorded 
during experiments with BMIs and demonstrate that the average 
timescale agrees well with the experimentally selected bin-width, 
though the timescales of individual units vary. The effect of 
the refractory period is assessed in Section Refractory Effects. 
Most importantly. Section BMI and Training Effects demon- 
strates that the timescales did not change significantly when 
switching to brain control or with training. During locomotion, 
the analysis identified the units whose timescales varied with 
the speed of walking — in agreement with the spectral analysis 
(Section Timescales During Walking). We conclude in Section 
Conclusions and Discussion, by discussing the proposed signal- 
independent timescale analysis (SITA) and its significance for 
investigating the timescales of neural-rate coding under different 
task conditions, and how they are related to the spectral content 
of the encoded signals and the refractory properties of the unit. 

MATERIALS AND METHODS 
MODELING AND ANALYSIS METHODS 
Rate-modulations Signal-to-noise ratio of DSPP 

The Poisson process is the simplest point process in which the 
probability of an event — a spike in the case of neural activ- 
ity, is assumed to be independent of the history of the spike 
train. Inhomogeneous Poisson processes are the simplest point 
processes that can describe rate modulations, i.e., the probabil- 
ity of a spike may change with time (Snyder, 1975; Cox and 
Isham, 1980; Johnson, 1996; Gabbiani and Koch, 1998). Here we 
consider DSPP, which is a subclass of inhomogeneous Poisson 



processes whose rate is modulated by stochastic signals (Snyder, 
1975; Cox and Isham, 1980; Johnson, 1996; Gabbiani and Koch, 
1998). Thus, the statistics of spike trains generated by DSPPs are 
determined by two factors: stochastic changes in instantaneous 
spike-rate and Poisson probability of spike occurrence. 

Considering long intervals of movements, the stochastic sig- 
nals that modulate the instantaneous rate, and thus the instan- 
taneous rate itself, are assumed to be stationary. This is justified 
since we consider long sequences of movements (reaching or 
stepping) that are not synchronized to any event, and thus may 
start at an arbitrary phase of the movement. This approach is 
illustrated in Figure 1 and contrasted with synchronized aver- 
aging, using a simple example where the rate is assumed to 
encode both the position (along 1 -dimension) during a reach- 
ing movement and an independent hidden signal. Figure lA 
illustrates the averaging method, where the same reaching move- 
ment (upper panel) is repeated, and the resulting spike-rates 
(lower panel, depicting a sample of 5 realizations) are synchro- 
nized to movement initiation. The resulting synchronized rate is 
a non-stationary process with time-varying mean (solid black line, 
averaged over a simulated sample of 500 realizations). Figure IB 
illustrates the proposed approach, showing a sample of five move- 
ment sequences (upper panel) starting at random phases (for 
illustration, each sequence includes four consecutive reaching 
movements between random targets). The resulting spike-rates 
(lower panel) are not synchronized and thus generate a station- 
ary process with constant mean and variance (solid and dashed 
black lines, respectively, based on a sample of 500 realizations). 
The rate-variance (Figure IB, lower panel) captures the contri- 
butions of both the movement (whose variance is depicted as 
black dashed line in Figure IB, upper panel), and the hidden sig- 
nal (whose variance is depicted as black dashed line in Figure lA, 
lower panel). 

Thus, the instantaneous spike-rate k{t) = -|- X(f) is 
assumed to be a stationary stochastic process with mean rate 
Ao and zero-mean rate modulations k(t). In realizations of 
DSPP, the distribution of spike-counts, Nj, in bins of size T is 
determined by the integrated spike-rate during the bin: 

Ar(f)= / X{a)da, (1) 

Jt-T/l 

and its statistics are given by (Snyder, 1975; Zacksenhouse et al., 
2007): 

E[NT^ = E[At] 
VarlNr] = VarlAr] + E[At] (2) 

The last equation can be interpreted as a decomposition of the 
total variance in the spike-counts into the variance of the mod- 
ulated rate, VarlAr], and the variance of the noise -E[Aj] = 
E[Nt] (Zacksenhouse et al, 2007). A similar decomposition was 
derived for general point processes from the total variance equa- 
tion (Churchland et al, 2011), however the restriction to DSPPs 
assures that the meaning of the two sources of the neural vari- 
ance is clearly defined (as further detailed in Section Conclusions 
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and Discussion), and facilitates explicit analysis of the effect of the 
refractory period, as detailed in Section Dead-time Effect. 

Considering the modulated rate as the signal, the RM-SNR is 
defined as: 



SNRt 



VariAr] 
E[At] ' 



(3) 



Using Equation (2), RM-SNR can be estimated from the mean 
and variance of the binned spike-counts as: 



VarlNr] - E[Nt] ^ ^ 

SNRnt = = _F — 1 

^ E[Nt] 



(4) 



vihere the hat denotes estimation, and _F is the Fano factor, defined 
as the ratio between the variance and the mean of the spike-counts 
(Dayan and Abbott, 2001). Thus, for the simple DSPP case, the 
KM-SNR is the deviation of the Fano-factor from its value for the 
homogeneous Poisson process, where _F = 1. Positive RM-SNR 
reflects rate modulations and characterizes irregular spike trains 
(i.e., spike trains whose variance is larger than their mean and 
thus larger than the variance of a homogeneous Poisson process 
with the same mean). Expressing the SNR is advantageous for fur- 
ther analysis of the timescales, as detailed in Section Timescales 
Analysis. 

Bandwidth effect 

The relationship between RM-SNR and the spectrum (co) of the 
instantaneous rate X(f) is analyzed in Appendix A. Equation (Al) 
shows that: 



SNRt - 



1 1 r 
XoTljt 



(5) 



where |Gw7 (ft')| = T Fourier transform of the 

rectangular window of duration T (Bendat and Peirsol, 2000). 
Furthermore, under the assumption that the spectrum of the 
instantaneous spike-rate S^(tt)) is band limited in the range 
[— (Wmax, Winax]> Appendix A derives two asymptotes for BM-SNR 
that hold at short and long bin-widths, respectively. For short 
bin-widths T, i.e., when (Wmax ^ In /T, Equation (A2) indicates 
that RM-SNR increases linearly with the bin-width. For long bin- 
widths r, i.e., when Wmax 5?> 2:^/7, Equation (A3) indicates that 
RM-SNR saturates. 

Dead-time effect 

During absolute refractory period, or dead-time r^, the instanta- 
neous firing rate is zero regardless of the modulating signals. The 
effect of the dead-time on the estimated RM-SNR is detailed in 
Appendix B, and, assuming that the dead-time is short compared 
to the dynamics of the instantaneous rate, an exact expression for 
SNRnj is derived (Equation B4). 

To gain more insight into the effect of the dead-time, the 
estimated SNRmj- is approximated (Appendix B, Equation B9) by: 



SNRNMd) 



SNRt 

[1 +TdAo] 



rgXp (2 -I- TdXp) 



(6) 



The approximation indicates that SNR^.,. estimated from a real- 
ization of a dead-time modified DSPP is a scaled and downward 
shifted version of SNRt associated with the original, dead-time 
free, rate. The scaling reflects the effect of the dead-time on the 
integrated rate (Equation B6), and the shift reflects the additional 
effect on the variance of the spike-counts (Equation B3). 

Simulations 

The analysis is demonstrated on simulated realizations of point 
processes with known characteristics. Specifically, the instanta- 
neous rate was generated as a stationary stochastic process by: 
(a) passing a unit variance white Gaussian noise through a 2nd 
order Butterworth filter with known cut-off frequency fcut, (b) 
scaling the resulting signal to achieve the desired variance of 
the instantaneous rate, and (c) adding a constant mean rate Xq. 
Finally, simulated spike trains were generated from the resulting 
instantaneous rate using the inverse distribution function tech- 
nique (Johnson, 1996). Spike trains were generated as realizations 
of DSPPs, doubly stochastic Gamma processes (DSGP; Fujiwara 
et al, 2009) with constant shape parameter k, or dead-time 
modified DSPPs or DSGPs with fixed dead-time. In summary, 
the simulated spike trains were characterized by five parameters: 
bandwidth defined by the cut-off frequency fcut of the low-pass 
filter, mean (Xq) and variance (Vx) of the instantaneous rate, 
dead-time rj, and the shape parameter for DSGP (/c, /<" = 1 for 
DSPP). 

We note that for homogeneous Gamma processes with k > I, 
the probability of firing after a spike is reduced, resulting in rel- 
ative refractory period, and the variance of the spike-counts is 
lower than for the homogeneous Poisson process with the same 
mean. In contrast, for homogeneous Gamma processes with k < 
1, the probability of firing after a spike is enhanced, and the 
variance of the spike counts is larger than for the homogeneous 
Poisson process with the same mean. 

Timescales analysis 

While SNRt, and hence SNRnj, increase with the bin-width, the 
update rate decreases (Wu et al, 2006) and the reaction time 
increases (Cohen and Newsome, 2009). We quantify the trade-off 
between SNRmj- and update-rate by evaluating their product, i.e., 
the ratio SNRnj-ZT, as a function of the bin- width T, and charac- 
terize the timescale of neural encoding as the bin-width at which 
this ratio saturates or peaks. This is demonstrated in Figure 2, 
which depicts the variations of SNRmj-/T with the bin- width 
for simulated realizations of DSPPs and DSGPs generated from 
the same realization of the instantaneous rate. The instantaneous 
rate was generated with mean Xq = 15 s^\ variance Vx = 48 s^^, 
and cut-off frequency /c„( = 1.0 Hz. The realizations analyzed in 
Figure 2A are dead-time free, while the realizations analyzed in 
Figure 2B include dead-time = 1 ms. 

Timescale analysis of DSPPs. The estimated SNRn^/T from both 
the dead-time free DSPP (Figure 2A, blue line) and dead-time 
modified DSPP (Figure 2B, blue line) closely match the val- 
ues computed from the rate (Equations 3 and B4, respectively, 
black lines). The curves depict the expected effects of the spec- 
tral content and dead-time. At long bin-widths, SNRt should 
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FIGURE 2 I Timescale analysis of simulated realizations of doubly 
stochastic Poisson and Gamma processes (DSPP, DSGP) (A) and 
dead-time modified {tj = 1) DSPP and DSGP (B). DSGPs are 
simulated with /c = 0.85, k = ^ (i.e., DSPP), and /c = 1.15. All 
simulations were generated from the same instantaneous rate function 
generated with mean >.o = 15s"\ variance Vi = 48s"^, and cutoff 



frequency fcut= 1-0 Hz. Estimated SNR^j/T (Equation 4) are compared 
with the values computed from the rate lusing Equation (3) for the 
dead-time free DSPP in (A), and Equation (B4) for the dead-time 
modified DSPP in (B)]. Analysis is conducted at bin-widths of 30, 40, 
50, 60, 70, 80, 90, 100, 110, 120, 130, 140, 150, 175, 200, 250, 300, 
500, 750, and 1000 ms. 



reach a constant asymptote (Equation A3), so SNRj/T should 
indeed decrease with the bin-width. The additive bias in estimat- 
ing SNRnj- from dead-time modified DSPP is constant (Equation 
6), independent of the bin-width, so, even when the spike-trains 
include a refractory period, the effect on the estimated SNRmj-/T 
would diminish at long bin-widths, as is evident in Figure 2B. 
At the other extreme, i.e., when the bin-widths are short enough 
for the linear asymptote (Equation A2) to hold, SNRj/T should 
be approximately constant, in agreement with the estimated 
SNRmt/T for DSPP (blue curve. Figure 2A). Hence, the short 
bin-widths at which SNRnj/T saturates are indicative of the 
timescales above which the slower than linear increase in SNR 
does not justify the reduction in update rate. 

Dead-time introduces a constant bias in estimating SNRn-^- 
(Equation 6), which becomes increasingly significant as SNRj 
diminishes at short bin-widths. Since the bias is negative, the 
estimated SNRn-p/T becomes smaller as the bin-width becomes 
shorter, in agreement with the short bin-width trend of the blue 
curve in Figure 2B. The effect of the dead-time at short bin- 
widths limits the benefit of reducing the bin-width, and results in 
a peak in the SNRm-i-/T curve. The bin-width at which this peak 
occurs maximizes the trade-off between SNRm-j- and update rate. 
We also note that below this bin-width the distortion due to the 
dead-time increases significantly. 

Timescale analysis of DSGPs. Applied to doubly stochas- 
tic Gamma processes (DSGPs) with scale parameter (k) 
larger/smaller than 1, the estimated SNRm-^-ZT under/over 



estimates the value computed from the rate, as demonstrated 
in Figure 2. In realizations of Gamma processes with k > I, the 
occurrence of a spike has an inhibitory effect, resulting in rela- 
tive refractory period. The estimated SNRn-i-/T curve in such a 
case (red line. Figure 2A) or when an absolute refractory period is 
also included (red line. Figure 2B) depict a peak similar to that for 
dead-time modified DSPP. As noted above (Section Simulations) 
realizations of homogeneous Gamma processes with k > I are 
characterized by spike-count variance that is smaller than the 
mean. The positive estimated SNRnj- (which implies variance 
larger than mean, see Equation 4) reflects the contribution of the 
variance of the stochastic rate. 

In realizations of DSGPs with k < I, the occurrence of a spike 
enhances, rather than inhibits, the probability of firing. The esti- 
mated SNRnt/T curve in such a case (green line. Figure 2A) 
increases sharply as the bin-width become shorter without sat- 
urating. For the specific parameters selected here, the excitatory 
recovery after a spike dominates the inhibitory effect due to the 
dead-time, and a similar increase is observed even with dead-time 
(Figure 2B). 

Timescales from SNRt^^/T curves. In summary, the analysis and 
simulations indicate that: (i) DSPPs are characterized by decreas- 
ing SNRmj/T curves that saturate at short bin-widths. In those 
cases, the bin-widths at which the SNRn-^/T curve saturates are 
indicative of the timescales above which the slower than linear 
increase in SNR does not justify the reduction in update rate, 
(ii) dead-time modified DSPP or DSGP with relative refractory 
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period (k > I) are characterized by SNRnj./T curves that peak 
at some bin- width. The bin-width at which the SNRmj-/T peaks 
maximizes the trade-off between SNRm^- and update rate, while 
limiting the distortion due to the refractory effect, (iii) DSPGs 
with excitatory recovery (k < I) are characterized by decreasing 
SNRnj/T curves that remain convex even at short bin- widths 
without saturating. In those cases, the appropriate time-scale 
is not evident from the proposed analysis, since the higher 
SNRn-j./T at short bin -widths reflect recovery effects rather than 
actual rate modulations. Further analysis is needed to select the 
bin-width that optimizes the trade-off between high SNRmj/T 
and small distortion due to the excitatory recovery effect. 

Variance, bandwidth and dead-time effects. The effects of the 
variance, bandwidth and dead-time are further evaluated in 
Figure 3, based on the approximation derived in Equation 6. 
Figure 3A establishes that the approximated SNRn^/T (dashed 
black line) captures well the main effect of the dead-time. Hence, 
the above effects are evaluated by approximating SNRn^ directly 
from the simulated stochastic rate functions. Unless otherwise 
specified, the nominal parameters of the stochastic rate functions 
are: mean rate Xq = 15s^\ variance Vi = 48 s^^, bandwidth 
= 1 Hz, and dead-time = 0.7 ms. 

Figure 3B demonstrates that the major effect of increasing 
the variance (from Vi = 24 s"^ to Vx = 48 s"^ and Vi = 72 s"^) 
is to increase the SNR. This is also accompanied by a shift in 
the peak location toward shorter bin-widths, due to the inter- 
action between SNRj and dead-time (Equation 6). Figure 3C 



demonstrates that as the cutoff frequency f^ut increases from 
0.5Hz to 1.0 Hz, and 1.5Hz, SNR^.j./T peaks at shorter bin- 
widths, capturing well the reduction in the relevant timescales. 
Figure 3D demonstrates that as the dead-time increases to 0.7 
and 1.4 ms, SN_R]Vj /r peaks at longer bin- widths. This trend 
reflects the increasing effect of the dead-time on the magnitude 
of the negative bias in Equation (6) (see note following Equation 
B9). Thus, for dead-time modified DSPP, the bin-width at which 
SNRn-i'/T peaks reflects the conflicting effects of the bandwidth 
of the encoded signals and the refractory period. 

EXPERIMENTAL METHODS 

The analysis tools developed above to estimate the KM-SNR of 
spike-trains and evaluate their timescales are applied to spike- 
trains recorded form cortical units in two experiments with 
Monkeys: (a) BMI experiments and (b) Locomotion experiments. 

Ethics statement 

All experiments were conducted with approved protocols 
from the Duke University Institutional Animal Care and Use 
Committee and were in accordance with the NRC/NIH guidelines 
for the Care and Use of Laboratory Animals. 

BMI experiments 

As detailed in Carmena et al. (2003), the BMI experiments were 
conducted with macaque monkeys and included three control 
modes: (i) pole control, (ii) brain control with hand move- 
ments (BCWH), and (iii) brain control without hand move- 
ments (BCWOH). During pole control the monkey controlled the 





FIGURE 3 I Approximated SNRi\jj/T for dead-time modified DSPP 
(Equation 6) in comparison with the S/VRju^/^ computed from the 
rate (Equation B4) or estimated from the spilce-counts (A). The 

approximation is used to demonstrate tine effects of the variance 



(B) bandwidth (C) and dead-time (D). Unless otherwise specified, the 
nominal parameters of the stochastic rate functions are: mean rate 
Xo = 15s"\ variance Vj. = 48 s"^, bandwidth fcut = 1 Hz, and dead-time 
= 0.7 ms. 
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cursor on the screen using a hand held pole. Data from pole con- 
trol was used to train a linear filter to predict the velocity from the 
neural activity. During brain control the cursor was controlled 
to follow the velocity predicted from the recorded neural activ- 
ity using the previously trained linear filter. Initially, the monkey 
continued to move its hand during brain control (BCWH), but 
starting at the 7th session, the monkey stopped moving the hand 
(BCWOH). Spike trains were recorded from single and multi- 
units (see Section Refractory Effects for more details) in the 
primary motor area (Ml), the pre-motor area (PMd), the supple- 
mentary motor area (SMA), and somato-sensory area (SI). Here 
we focus on significantly modulated units (n = 163), i.e., those 
units for which the null hypothesis that their ISI distribution was 
generated from a homogeneous Poisson process can be rejected at 
a = 0.05 significance level (Zacksenhouse et al., 2007). 

Locomotion experiments 

During the locomotion experiments, rhesus macaques either 
stood or walked bipedally on a treadmill in three different speeds, 
12.5, 25, and 50cm/s, as detailed in Fitzsimmons et al. (2009). 
Spike trains were recorded from single and multi-units (66 and 
34%, respectively) in regions associated with the representation 
of the lower limbs in the primary motor (Ml ) and somatosensory 
(SI) cortices. 



RESULTS 

SNR DURING BMI EXPERIMENTS 

Figures 4A,C depict the mean SNRn-,. and mean normalized 
SNRnj. estimated from spike-trains of significantly modulated 
cortical units (« = 163) recorded during a BMI experiment, as 
a function of the bin-width T. Assuming ergodicity. Equation 
(4) was used to estimate the SNRn^^, of each unit from its 
spike-counts statistics in non-overlapping windows of 2 min. 
The mean SN_RjVj- in Figure 4A was derived by first computing 
ensemble mean (over n = 163 units) at each window, and then 
computing the mean and standard deviation over 10 consecu- 
tive windows in each of the experimental modes (pole control 
and brain control with and without hand movements). Thus, 
error bars in Figure 4A depict the standard deviations of the 
ensemble-means over the 10 windows of time. Variations across 
the ensemble of units are described in Figure 4C, which depicts 
the ensemble mean and standard deviations of the normalized 
SNRnj - The normalized SNRn.^. for each unit was computed 
by first averaging across all windows, and then normalizing by 
a scaling factor. The scaling factor was selected so the nor- 
malized SNRmj- at bin-width of 300 ms is the same for all 
units and equals the original ensemble mean. Thus, error bars 
in Figure 4C depict the standard deviations of the normalized 
SNRmj across the ensemble of significantly modulated units 
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FIGURE 4 I Mean estimated SNRrj^ (A) and normalized S/VR/v^ (C), 
and the corresponding S/Vf?/v^/r (B,D) of neural activity recorded 
from significantly modulated units (n= 163) during BIVII experiments. 

BMI experiments included pole control, brain control with hand 
movements (WH), and without hand movements (WOH). SA/R/v^ was 
computed from the spike-train of each unit in non-overlapping windows of 
2-min. (A,B) are based on first averaging over the units at each window, 
and then computing the mean and standard deviation over 10 consecutive 
windows in each of the experimental mode. Hence, error-bars in 



(A,B) mark the standard deviations across 10 consecutive windows. (CD) 
are based on first averaging the SNRfjj of each unit across all windows, 
and then scaling it so the normalized SNRfj^ at bin-width of 300 ms is the 
same for all units and equal the original ensemble mean. Ensemble mean 
and standard deviation of the normalized SNRnj were computed across 
units. Hence, error bars in (CD) mark the standard deviations across the 
ensemble of units (resulting in zero standard deviation at 300 ms, which 
was used for normalizations). Dashed lines in (A) depict the linear fits for 
short bin-widths, over the range 50-1 50 ms. 
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(resulting in zero standard deviation at 300 ms, which was used 
for normalizations). 

Figures 4A,C demonstrate that the mean estimated SNRnj- 
increases approximately linearly at short bin-widths and 
approaches a constant level at long bin-widths, in agreement with 
the expected asymptotes (Equations A2 and A3). Another strik- 
ing feature is that at each bin-width, the mean SNRmj- in brain 
control is always higher than that in pole control; and is high- 
est in brain control without hand movements. This extends our 
previous results (Zacksenhouse et al., 2007) for T = 100 ms — the 
bin-width used for operating the BMI. 

TIMESCALE DURING BMI EXPERIMENTS 

Figures 4B,D demonstrate that the mean SNRn-p/T and mean 
normalized SNRn-^./T increase at short bin- widths, reach a wide 
peak around 100 ms, and decrease at long bin-widths. This pat- 
tern agrees well with the analysis and approximation of the 
SNRnj./T of dead-time modified DSPP, and with estimations 
derived from the simulated spike-trains of dead-time modi- 
fied DSPPs (or DSGPs with k < I, with or without dead-time) 
depicted in Figure 2. 

Focusing on individual units from PMd and Ml, the SNRn-^ /T 
curves shown in Figure 5 demonstrate two typical patterns: (i) 
curves with a wide peak around 100 ms (Figures 5 C,D), simi- 
lar to the average curves in Figures 4B,D and to the curves for 
dead-time modified DSPP or DSGP with k < 1, in Figure 2, or 
(ii) decreasing curves (Figures 5A,B), similar to the curves for 
DSPP or DSGP with k > I in Figure 2. The second type could 
be further divided depending on whether the curve saturates at 
short bin-widths as for DSPPs, or remains convex as for DSGPs 
with K > I. 

Units were classified as exhibiting a wide peak around 100 ms 
if the peak SNRnj-/T was in the range of 50-150 ms, and as 



exhibiting decreasing SNR^j./T if the peak was below 50 ms. 
Note that SNRmj-/T of few units peaked at bin-widths longer 
than 150 ms, and were not included in either group. To assure 
proper peak-detection, we restricted this classification and fur- 
ther analysis of the timescales to the significantly modulated units 
whose maximum SNRnj/T during pole control was larger than 
0.5 s"^ (« = 109). Of those units, 59.7, 64.2, and 61.5% exhib- 
ited a wide peak around 100 ms during pole control, BCWH 
and BCWOH, respectively, while 25.7, 21.1, and 22.0% exhib- 
ited decreasing curves, respectively. The second group was further 
divided depending on whether the decreasing curve was convex at 
short bin-widths, i.e., whether M2 < (ui -|- M3) /2whereMi, U2, M3 
are the mean SNRn-^./T in the bin-width intervals 30-50, 60-80, 
and 90-100 ms, respectively. About 15% of the curves were con- 
vex at short bin-widths (17.4, 11.0, and 15.6%, in pole control 
BCWH and BCWOH, respectively). In those cases, the time-scale 
analysis is not conclusive, and is either not appropriate since the 
recovery function is dominated by excitatory effects (as clarified 
in Section Timescales Analysis) or should be extended to shorter 
bin-widths. 

The SNRm'i-/T estimated form spike-trains recorded during 
BMI experiments (Figures 4B,D, 5) vary with the bin-width in a 
similar way to the SNRn-p/T of simulated spike-trains (Figure 2). 
Decreasing curves, which saturate at short bin-widths (Figure SB, 
pole and BCWH) agree well with the assumption that the spike- 
trains are realizations of DSPP while decreasing curves that 
remain convex at short bin-widths (Figure 5A, pole and BCWH) 
agree well with the assumption that the spike-trains are real- 
izations of DSGP with K < I. Curves that peak (Figures 5C,D) 
agree well with the assumption that the spike trains are realiza- 
tions of dead-time modified DSPP or DSGP with /c > 1 (with 
or without dead-time), which are characterized by relative, rather 
than absolute, refractory period. 




FIGURE 5 I Estimated SNR^j/T as a function of the bin-width for four representative cortical units in PIVId (A,C) and Ml (B,D), during pole control, 
and brain control with and without hand movements, WH and WOH, respectively. 
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Refractory effects 

To assess the potential contribution of refractory effects, we 
estimated the refractory interval from the ISI distribution. The 
minimum ISI in all units was 0.675 ms, and most units (95% 
of the significantly modulated units) had minimum ISI of less 
than 1 ms. Since a single ISI is prone to measurement errors, we 
estimated the dead-time as the interval which was exceeded by 
99% of the ISIs in at least one mode of operation. Of the 109 
units described above, 50% had dead- time greater than 1.6 ms, 
and thus could be considered single-units (Hatsopoulos et al., 
2004). The percent of units having decreasing curves was sim- 
ilar, independent of whether the dead-time was short or long 
(below or above 1.6 ms); but the percent of units having time- 
scales longer than 150 ms was larger for units with long dead-time 
(22.4% on average, across all modes, compared with 8% for 
units with short refractory interval). Furthermore, the bin- widths 
at which SNRmj/T peaked had a moderate positive correlation 
with the product of the dead-time and the mean rate (coeffi- 
cient of correlation of 0.4, 0.42, 0.48 in pole control BCWH and 
BCWOH, respectively). Both the positive correlation and its mod- 
erate value are in agreement with Equation (6), which indicates 
that this product reduces the SNRn^. (and hence shifts the peak of 
the SNRn-i'/T toward longer bin-widths, as clarified in Section 
Timescales Analysis), but that SNRn,. is also affected by other 
variables (i.e., SNRj). 

BMI AND TRAINING EFFECTS 

Figures indicates that the bin- width at which the SNRm-^-ZT 
peaks remains approximately the same even after switching to 
brain control. Specifically, the SNRn-i^/T curves in pole control, 
BCWH and BCWOH peak at: 30, 40, and 30 ms for the PMd unit 
in panel (A), 130, 150, and 140 ms for the PMd unit in panel (B); 
30, 40, and 30 ms, for the Ml unit in panel (C); and 150, 150, 
and 110 ms for the Ml unit in panel (D). The distribution of 



the bin- widths at which the estimated SNRn.j./T peak and how 
they change when switching to brain control is summarized in 
Figure 6. The bin- width at which the estimated SNRnj-/T peaks 
during pole control is similar to (within 50 ms of) the bin-width 
at which it peaks during brain control for most of the units (80.7 
and 67.9% for BCWH and BCWOH, respectively). 

The spectral content of the spike-counts (in 100 ms bins) of 
the four units whose timescales were analyzed in Figure 5, are 
depicted in Figure 7. Specifically, the power spectrum density 
(PSD) was computed from the spike-counts in bins of 100 ms 
after subtracting the mean and expressed in dB. The main effect 
of the transition to brain control is the increase in power, which is 
higher in brain control than in pole control and is highest in brain 
control without hand movements. This is consistent with the 
higher RM-SM? in brain control compared to pole control. The 
overall bandwidth of the activity of each unit remains relatively 
similar across the different BMI control modes [even when, as in 
panel (D), the PSD of the activity during pole control portrays a 
small peak which is surpassed by the higher overall power during 
brain control] — consistent with the relatively similar shapes and 
peak locations of the different SNRisij./T curves for each unit in 
Figure 5. 

Figure 8 demonstrates that the shape of the SNRnj-/T curves 
is also invariant to training. At later sessions, as training pro- 
gresses, the SNR at each control mode and bin-width decreases (as 
detailed in Zacksenhouse et al., 2007) for bin-width of 100 ms). 
However, the effect of bin- width is similar and the SNRn;./T 
curves peak at a similar bin-width in all sessions. 

TIMESCALES DURING WALKING 

Figures 9, 10 depict the estimated SNRn,j /T and the spectrum of 
representative cortical units in the motor leg area while the mon- 
key was standing (black) or walking at increasing speeds (12.5, 
25, and 50 cm/s). The SNRmj./T curves in Figures 9A,B peak at 
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FIGURE 6 1 Distribution of timescales during BMI experiments. 


control without hand movements (WOH, B). Only significantly modulated 
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FIGURE 7 I Power spectrum density (PSD) of the neural activity analyzed in Figure 4. Tlie spectrums were computed from spike-counts in 100 ms bins, 
witli frequency resolution of 0.039 Hz. 
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FIGURE 8 I Mean estimated SA/R/v^/T across the population of 
significantly modulated cortical units recorded in three different BMI 
sessions during pole control (PC), and brain control with and without 
hand movements (BCWH and BCWOH, repectively) as a function of 
the bin-width. Tine three sessions include: an early session (3rd day); a 
mid session (7th day), and a late session (10th day). The 7th day is the first 
one in which BCWOH was performed, and in subsequent sessions BCWH 
was not perfomrned. 



progressively shorter bin-widths as the speed of walking increases. 
During standing, the peak SNRnj-/T is shallower and appears at 
longer bin-widths than during walking. Figures 10A,B indicate 
that the spectrums of the spike-counts (in 100 ms bins) of these 
spike-trains are dominated by single spectral lines that appear 



at higher frequencies as the speed of walking increases. During 
standing, the dominating spectral line appears at very low fre- 
quencies and is smaller compared to the spectral lines during 
walking. Thus, this simple case, in which the spectrums are dom- 
inated by single spectral lines, helps to demonstrate the expected 
effect of changes in the spectral content of the modulating sig- 
nals on the timescale. However, the timescale emerges from the 
interaction between the spectral content and the refractory prop- 
erties of the individual unit. Hence, even though the dominating 
spectral lines appear in the same frequencies (Figures 10A,B) the 
timescales are unit specific (Figures 9A,B). 

In contrast, the neural activity from other units in the leg area 
are characterized by SNRnt/T curves that either decrease with- 
out peaking (Figure 9D) or depict a peak that does not shift as 
expected with the speed of walking (Figure 9C). In both cases 
the neural activity during standing is characterized by higher 
SNRmj-/T than during walking. The corresponding spectrums 
(Figures 10C,D) indicate that the power of the spike-counts 
recorded from these two units is indeed larger during standing 
than during walking. The spectral lines associated with the speed 
of walking (at the frequencies in which they appear in the spec- 
trums depicted in panels A and B) are very small. Finally, for 
both units the neural activity recorded during running has more 
power in the low- frequency range than during walking. This may 
explain the unexpected shift in the peak of the SNRj/T toward 
longer rather than shorter bin-widths when the Monkey is run- 
ning compared to walking (Figure 9C). The convex shape of the 
SNRn-i'/T curve in Figure 9D even at short bin- widths suggests 
that the activity of this unit is dominated by excitatory recovery 
effects, and the timescale is not conclusive. Nevertheless, since the 
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FIGURE 10 I Power Spectrum Density (PSD) of the neural activity analyzed in Figure 9. PSDs were computed from the binned neural activity with 100 ms 
bins and frequency resolution of 0.039 Hz. 



SNRn.j~/T curves are similar during walking and running (inde- 
pendent of the speed) but much higher during standing, it can 
be concluded that signals associated with the frequency of walk- 
ing have little effect on the spike-rate of these two units, and that 
their activity is more strongly modulated during standing. 



Additional analysis (not shown) indicates that the mean 
SNRnt/T and mean normalized SffRMj-/T across significantly 
modulated units with peak above 0.5 s^' (n = 48) depict a wide 
peak around 50-1 10 ms, in most cases. The only exception is the 
normalized SNRn-^-ZT during standing, which is highest at the 
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shortest bin-width analyzed (25 ms). Around the peak, the mag- 
nitude of the SNRm-i-/T are similar independent of the speed so 
the curves overlap each other. At longer bin-widths SNRnj/T 
during slow walking is higher than during fast walking or run- 
ning, and is highest during standing. 

CONCLUSIONS AND DISCUSSION 
SNR ESTIMATION 

We suggest a computational framework for estimating the SNR in 
spike trains without relying on any assumption about the encoded 
signals. The suggested approach captures all the encoded sig- 
nals, including hidden and internal signals, and not only those 
controlled experimentally. The proposed approach has two lim- 
itations: (i) it accounts also for possible task irrelevant rate 
modulations, and (ii) the estimation method is derived under the 
assumption that the spike-trains are realization of DSPP. The first 
limitation is a necessary side-effect of capturing the hidden or 
internal signals. The second limitation is partially relaxed by con- 
sidering dead- time modified DSPP and doubly stochastic Gamma 
processes. 

The dependency of RM-SM^ on bin-width was first analyzed 
for DSPP and two asymptotes were derived for bin-widths that are 
either very short or very long compared with the inverse of the 
bandwidth. In particular, it was shown that RM-SNR increases 
linearity at short bin-widths and reaches saturation at long bin- 
widths. Extending the analysis to dead-time modified DSPP, it was 
shown that the estimated KM-SNR is reduced by both scaling and 
shifting factors that do not depend on the bin-width. The bias 
effect becomes more significant as the RM-SNR becomes smaller 
at shorter bin-widths. Using simulations of doubly-stochastic 
Gamma processes, we also evaluate other recovery effects in which 
the occurrence of a spike temporarily and partially inhibits (as in 
relative refractory period) or enhances the probability of firing. 

INHOMOGENEOUS vs. DOUBLY STOCHASTIC POINT PROCESSES 

Doubly stochastic point processes are inhomogeneous point pro- 
cesses with stochastic rate function. The alternative model for 
describing rate modulations is the inhomogeneous point pro- 
cess with a deterministic rate function (Cox and Isham, 1980; 
Kass and Ventura, 2001; Nawrot et al, 2008). The determin- 
istic model is based on the assumption that the rate function 
depends solely on the experimentally controlled conditions, and 
thus that it can be estimated from the post-stimulus time his- 
togram. This assumption might be valid for peripheral neurons 
that are affected almost exclusively by the experimentally con- 
trolled stimulus, but less for cortical neurons which may encode 
additional hidden signals, not directly under experimental con- 
trol. These signals may vary from trial-to-trial and thus can only 
be modeled as stochastic. During movements, motor planning 
(Churchland et al., 2006a,b) and prediction errors (Wolpert and 
Ghahramani, 2000; Shadmehr and Krakauer, 2008), for example, 
may vary from trial-to-trial, while in decision making the accu- 
mulated evidence (while potentially sharing a common trend), 
may vary from trial-to-trial (Churchland et al., 2011). In these 
cases, doubly stochastic point processes are essential for captur- 
ing the variance of all the signals that might be encoded in the 
neural activity. 



VARIANCE DECOMPOSITION 

Under the assumption that spike-trains are realizations of DSPPs, 
the variance of the spike-counts is decomposed into two terms 
associated with the signal and noise. A similar decomposition was 
developed in Churchland et al. (2011) where the resulting two 
terms were described as the variance of the underlying "inten- 
sity command" and the variance of the neural activity given 
that command. That decomposition was applied to uncover how 
trial-to-trial neural variability changes during decision making 
and reveal the underlying nature of neural processing. WhUe the 
application and some of the analysis details differ (as discussed 
below), we share the most important hypotheses: (i) the neu- 
ral activity may represent hidden signals, i.e., signals that are 
not directly controlled and hence may vary from trial to trial, 
(ii) decomposition of the variance of neural activity facilitates 
estimating the variance of the underlying rate, including modula- 
tions by those hidden signals, and (iii) variations in the variance 
of the underlying rate with time (or task) may reveal important 
information about the neural processing. 

The decomposition developed in Churchland et al. (2011) is 
based on the total variance equation, which holds for any dou- 
bly stochastic point process. The total variance is decomposed 
into the variance of the conditional expectation (the variance of 
the expected spike-count given the rate parameter and the point 
process variance (the expectation of the conditional variance). 
Under the restriction to DSPP made here, the first term equals 
the variance of the rate parameter, which is assumed to encode the 
modulating signals. Otherwise, the relationship of the first term 
to the variance of the rate parameter is more complex. In par- 
ticular, for dead-time modified DSPPs and for DSGPs, the first 
term is proportional to the variance of the rate-parameter but 
not equal to it (see Equations Bl and B6 for the dead-time mod- 
ified DSPP). Hence, the restriction to DSPP made here assures 
that the meaning of the first term is clearly defined. Furthermore, 
it enables explicit analysis of the effect of the refractory period 
on both terms. Another major difference is that the analysis in 
Churchland et al. (201 1 ) was applied to estimate the synchronized 
trial-to-trial variance in the underlying rate-process (intensity 
command) at each time along the trial. In contrast, we applied the 
analysis to estimate asynchronous variance in the underlying rate 
process. 

TIMESCALES OF NEURAL ACTIVITY 

As bin-width increases, SNR increases but update-rate is com- 
promised (Wu et al., 2006). This trade-off is captured by the 
ratio of SNR to bin-width, SNRj/T. Thus, the bin-width at 
which SNRj/T peaks optimizes the trade-off between SNR and 
update rate and is suggested as a criterion for characterizing the 
timescales of neural rate-coding. 

Our analysis suggests that the estimated SNRnj-/T is affected 
by both the spectral content of the encoded signals and the 
recovery properties of the unit. The effects of the bandwidth 
and absolute refractory period were analyzed theoretically and 
demonstrated via simulations of DSPPs and dead-time modi- 
fied DSPPs. Their conflicting effects give rise to a peak in the 
SNRnt/T curve, which characterizes the timescale of neural 
rate-coding. 
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DEVIATIONS FROM THE DSPP ASSUMPTION 

The DSPP assumption is violated when the firing rate depends on 
the history of the spike train. In these cases, the estimated SN_R]v.j- 
may differ from the actual SNRj, due to history-dependent mod- 
ulations of the instantaneous rate. We restrict the analysis to the 
simple, yet common case, of renewal point processes, in which 
the dependency on the history is captured by the time elapsed 
since the last spike. The effect of absolute refractory period was 
analyzed theoretically, and was shown to result in a peak in 
the SNRnj/T curve. Thus, the dead-time limits the benefit of 
reducing the bin-width, since the estimated SNRn-i- becomes 
increasingly distorted. The bin-width at which the peak occurs 
maximizes the trade-off between SNRmj- and update rate, and 
characterizes the timescale of neural rate-coding. 

Using simulations we also demonstrated that relative refrac- 
tory periods have similar effect, resulting in SNRnj-/T curves 
that are also characterized by a peak. So the timescales are sim- 
ilarly characterized by the bin-widths at which the SNRm-i-/T 
curves peak. In contrast, when the recovery period is charac- 
terized by enhanced, rather than inhibited, probability of firing, 
the SNRm-i-/T over-estimates (rather than under-estimates) the 
SNRt/T. Thus, as the bin-width decreases, SNRn.,-/T increases 
without peaking or saturating. In those cases, the appropriate 
time-scale is not conclusive from the proposed analysis, since 
the higher SNRmj-ZT at short bin-widths reflect recovery effects 
rather than actual signal modulations. Further analysis is needed 
to select the bin-width that optimizes the trade-off between high 
SNRn.,~/T and small distortion due to the excitatory recovery 
effect. 

TIMESCALES DURING BMI EXPERIMENTS 

Estimating RM-SNR from neural activity recorded during BMI 
experiments, we observe that while its magnitude depends on the 
BMI mode and on training, the timescale revealed by the peak of 
the SNRn-i-/T curve remains the same. 

The observation that KM-SNR increases when switching to 
brain control indicates that the variance of the rate increases — 
and thus that the variance of some of the encoded signals 
increases. The observation that the timescales remain invariant 
suggests that those signals are task- relevant, and we are currently 
investigating the hypothesis that they include state estimation 
and control signals, within the framework of optimal control. 
However, their exact nature and decoding remains an open issue. 

Most of the analyzed SNRn-i^/T curves (62% on average across 
the different control modes, n = 109) were characterized by a 
wide peak at 50-150 ms. For this class of curves (and for the 
15% of curves that peak at longer bin-widths), decreasing the bin- 
width below the peak may compromise decoding, unless the effect 
of the refractory period is accounted for. Indeed BMI applica- 
tions based on spike-counts used bin-widths in the range between 
30 and 100 ms (Wessberg et al, 2000; Serruya et al, 2002; Kim 
et al., 2008; Velliste et al., 2008), consistent with the range of 
timescales revealed here. In particular, the BMI experiments ana- 
lyzed here were conducted with 100 ms bins (Carmena et al., 
2003) within the wide peak of the estimated SNRmj-/T curve 
depicted in Figure 4B. Our analysis suggests that this selection 
provides a good trade-off between the SNR and update rate. 



Curves whose peak was below 50 ms were divided into those 
that were convex or concave at short bin-widths (15 and 8% 
on average, respectively). Concave shape at short bin-widths is 
expected from SNRn,^./T curves derived from DSPPs, where this 
ratio should saturate at short bin-widths as detailed in Section 
Timescales Analysis (and Figure 2A). In this case, the range of 
bin-widths at which the curve saturates is indicative of the rel- 
evant timescales for neural decoding. At longer bin-widths, the 
slower than linear increase in SNR does not justify the reduction 
in update rate. Convex SNRmj/T are expected when the unit has 
excitatory recovery period (in which the probability of firing is 
enhanced). As detailed above, the timescale analysis in this case is 
not conclusive. 

Since the relevant timescales vary across individual units, 
it might be beneficial to apply different bin-widths when bin- 
ning different units. This may explain the potential advantage 
of multi-resolution binning (Kim et al, 2005). Recent decod- 
ing methods are based on point process filters that operate at 
5 ms bins (Shanechi et al., 2013a,b). The BMI tested in Shanechi 
et al. (2013a) was based on multi-unit recording, where the dead- 
time effect might be negligible, and hence the timescales could 
be short. Furthermore, point process filters are based on estimat- 
ing the conditional probability of the spike-counts given the rate, 
which is computationally more efficient when the bin is small 
(and the number of spikes that may occur with significant prob- 
ability is limited). Hence, the choice of small bins may reflect 
computational constraints. The analysis conducted here suggests 
that unless recovery effects are negligible or accounted for, very 
short bin-widths may compromise the trade-off between SNR 
and update-rate. 

TIMESCALES DURING LOCOMOTION EXPERIMENTS 

The overall shape of SNRn;-/T curves estimated from spike trains 
recorded during walking is similar to those derived from the spike 
trains recorded during the BMI experiments. However, the peak 
location in the locomotion experiments varied with the task, at 
least for some units, and shifted to lower bin-widths when the 
speed of walking increased. This is in agreement with the expected 
change in the bandwidth of encoded signals associated with the 
experimentally controlled speed of walking. Indeed, units whose 
timescales vary in this way are characterized by PSDs dominated 
by a single peak at the frequency of walking. Thus, the locomotion 
test-case demonstrates the expected effect of changes in the spec- 
tral content of the modulating signals on the timescales. However, 
timescales depend not only on the spectral content but also on the 
refractory properties of the units. So even if the spectral content 
is similar the timescales are unit-specific. 

In summary, we developed a method for assessing the SNR and 
timescales of neuronal activity independent of the encoded sig- 
nals and observe that many units depict well-defined timescales 
of the order of tens to hundreds msec. The timescale emerges 
from the interaction of the recovery properties of the unit and 
the spectral content of the modulating signals. In particular, 
the variations of the estimated SNRmj/T at short bin-widths 
is related to the recovery properties of the unit. During a 
BMI experiment, the average timescale across the population 
agrees well with the experimentally selected bin-width, though 
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individual units depict different timescales. Furthermore, while 
SNR increases after switching to brain control and decreases 
with training, timescales remain similar. During locomotion, 
the analysis identified units whose timescales varied consistently 
with the experimentally controlled speed of walking, though 
the specific timescale reflected also the recovery properties of 
the unit. Hence, the proposed method provides a hypothe- 
sis free tool for investigating possible contributions of hidden 
signals, not under experimental control, to neural modula- 
tions and the effect of task conditions on their magnitude and 
timescales. 
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APPENDIX A: BANDWIDTH EFFECT 

The binning process is equivalent to passing the instantaneous 
spike-rate through a finite time integrator followed by down- 
sampling and spike-counting (see Bair et al., 1994 for a similar 
formalization of the recording of spike-trains with finite time res- 
olution). Specifically, the integrated spike-rate Aj(f) = XqT + 
Aj{t), where A j-(f) = Jl^jj^ X(cr)c?cr, is obtained by convolving 
the time-dependent component of the instantaneous spike-rate 
k{a) with a rectangular window of duration T: Ar(f) = X(f) (gi 
Wjit) [(g) denotes convolution, and Wj(f) is a rectangle signal 
of duration T (Bair et al., 1994)]. The effect of down-sampling 
will be considered later. For now we note that the spectrum of 
Aj(f), denoted by S^(tti), is related to the spectrum S-^(u)) of A.(f) 

by:S^(w)= \Gwt{m)^ Si{m),whtre\GwT{M)\ = T 
is the Fourier transform of the rectangular window (Bendat and 
Peirsol, 2000). The variance of A7(f) can be computed from 
the spectrum S^(ci)) as: VarlAjit)] = f^^S^{a>)da), while 
its mean is _E[A j(f)] = ElNjit)] = koT. Using Equation (3), the 
SNR of the binned spike-counts is given by: 



SNRi 



1 1 r 

XoTItt 



\Gv/-i;{u>)\ S-^{u>)du> (Al) 



Assuming that the spectrum of the instantaneous spike-rate 



Sj{w) is band limited in the range [- 



the SNR 



reaches two asymptotes for short and long bin-widths. For 
short bin-widths T, i.e., when comax <K In /T, the transfer 
function of the rectangular window |Gw'j (ft')| in the range 
[— Wmax, Wmax] Can be approximated by |Gw.^(0)| = T, so, from 
Equation (Al): 



SNRt ■■ 



1 



f 

J —C 



S-^(co)dco ■■ 



ko 



(A2) 



where Pz 



X ~ Jn f-'(T^ S-^((o)da) is the power of the time- 
dependent component of the instantaneous rate. Hence, at short 
bin-widths, the SNR of spike-counts generated from DSPP with 
band-limited instantaneous rate should increase linearly with the 
bin-width T. 

For long bin-widths T, i.e., Wmax 5S> Ztt/T, the Fourier trans- 
form of the rectangular window is concentrated around zero, 
so only low frequencies pass the averaging process and affect 
the average rate. Denoting by So the spectrum of the instanta- 
neous rate at low frequencies. Equation (Al) implies that at long 
bin-widths the SNR saturates at: 



SNRt 



1 Sg 



f 

J —I 



Gw.j (w)| dca ■■ 
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Binning is usually performed in non-overlapping windows, so 
the binned spike-counts are related to the samples of the inte- 
grated spike-rate taken at sampling interval T. At short bin- 
widths, the Nyquist frequency WNy = oisamplmg/'i- = n/T satisfies 
the Nyquist criterion WjVy > Mmax. Hence, there is no aliasing, 
and Equation (A2) holds, implying that the SNR of binned spike 
counts should increase linearly with the bin-width. At long bin- 
widths, the Nyquist criterion is not satisfied, and the side bands 



of the spectrum of the rectangular window | Gv\f.^.{w) | fold back 
into the main lobe. Nevertheless, as long as the spectrum is 
constant over most of the side-lobes, the total power in the 
main lobe (including all the folded back side-lobes) is the same 
as the integral computed in Equation (A3) (Bair et al., 1994). 
Hence the SNR of binned spike counts should saturates at long 
bin-widths. 

APPENDIX B: SNR OF DEAD-TIME MODIFIED DSPP 

Assuming the dead-time is short compared to the dynam- 
ics of the instantaneous rate, the statistics of the spike-counts 
depend on the dead-time modified integrated rate: Aj = 

ft-T/i i+^rj^o-) '^'^' ^^'^ given by [see Equations 3 and 9 in 
Vannucci and Teich (1981)]: 



E[Nt]=E[At] 

£[N^]=£[(Ar)']+£ ^ ^ ^ ^ ^^^^^^3 



'•t+T/2 



X{a) 



t-T/2 [l + rdX{(r)V 



da 
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(B2) 



So: 



Var[NT^ = Var[hT^+E 



Jt-T/2 ll + rdXia)r 



(B3) 



Substituting Equations (Bl) and (B3) in Equation (4), the esti- 
mated SNR from a dead-time modified DSPP is: 



Var 



SNRNMd) 



E[At] 



(B4) 



The SNR associated with the dead-time modified integrated rate, 
SNRt, can be defined in a similar way to the definition of SNRt 
for the original integrated rate (Equation 3): 



SNRt ■■■ 



Var [At] 
E[At] 



(B5) 



For rciX(a) <^ 1, the integral defining Aj can be approximated 

by: At = i-f-]^Xa ft-T/i Mcf)'icr = i^pi^' 5° statistics can be 
approximated by: 



E[At] 
Var[AT] 
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Hence, the SNRt can be related to SNRt by: 

SNRt 
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Using a similar approximation, the second term in Equations 
(B2-B4) can be approximated as: 
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SNRt 
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Substituting these approximations in Equation (B4) results in an 
approximated relationship between the estimated SNR from a 

dead-time modified DSPP and the SNR of a dead-time free DSPP Note that the second term increases with the dead-time both 



having the same instantaneous rate: 



absolutely and as a fraction of the first term. 
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